#include "KDM1DetectorConstruction.hh"
#include "KDM1Globals.hh"
#include "KDM1Detector.hh"
#include "G4Element.hh"
#include "G4Material.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4LogicalVolume.hh"
#include "G4ThreeVector.hh"
#include "G4PVPlacement.hh"
#include "G4SDManager.hh"
#include "G4UImanager.hh"
#include "G4ios.hh"
#include "G4NistManager.hh"
#include "G4MaterialPropertiesTable.hh"
#include <cmath>
#include <sstream>

KDM1DetectorConstruction::KDM1DetectorConstruction()
{
	;
}

KDM1DetectorConstruction::~KDM1DetectorConstruction(){;}

G4VPhysicalVolume *KDM1DetectorConstruction::Construct()
{
	// Define materials
	
	G4NistManager* nistMan = G4NistManager::Instance();

	G4Element* N = new G4Element("Nitrogen","N",7,14.01*g/mole);
	G4Element* O = new G4Element("Oxygen","O",8,16.00*g/mole);
	G4Material* Air = new G4Material("Air",15*g/6.022/pow(10,23)/cm3,2);
	Air->AddElement(N,70*perCent);
	Air->AddElement(O,30*perCent);

	G4Material* CsI = nistMan->FindOrBuildMaterial("G4_CESIUM_IODIDE");

	// Define world
	
	G4Box* worldSolid = new G4Box("worldSolid",2.5*m,2.5*m,5*m);
	G4LogicalVolume* worldLog = new G4LogicalVolume(worldSolid,Air,"worldLog");
	G4VPhysicalVolume* worldPhys = new G4PVPlacement(0,G4ThreeVector(0,0,0),worldLog,"world",0,0,0);

	// Define crystal block
	
	G4Tubs* crystalBlockSolid = new G4Tubs("crystalBlockSolid",0,0.95*m,25*cm,0,360*deg);
	G4LogicalVolume* crystalBlockLog = new G4LogicalVolume(crystalBlockSolid,Air,"crystalBlockLog");
	G4VPhysicalVolume* crystalBlockPhys = new G4PVPlacement(0,G4ThreeVector(0,0,25*cm),crystalBlockLog,"crystalBlock",worldLog,0,0);
	(void)crystalBlockPhys;

	// Define small crystals

	G4SDManager* SDMan = G4SDManager::GetSDMpointer();
	G4VSensitiveDetector* sensitiveDetector = new KDM1Detector("sensitiveDetector");
	SDMan->AddNewDetector(sensitiveDetector);

	G4Box* smallCrystalSolid = new G4Box("smallCrystalSolid",1.25*cm,1.25*cm,25*cm);
	G4LogicalVolume* smallCrystalLog = new G4LogicalVolume(smallCrystalSolid,CsI,"smallCrystalLog",0,sensitiveDetector);

	int crystalCount = 0;

	int smallCrystalDimension = 48;
	for(int i = 0;i < smallCrystalDimension;i++)
	{
		for(int j = 0;j < smallCrystalDimension;j++)
		{
			double xPos = (i-smallCrystalDimension/2+.5)*2.5*cm;
			double yPos = (j-smallCrystalDimension/2+.5)*2.5*cm;
			if(fabs(xPos) > 10*cm || fabs(yPos) > 10*cm)
			{
				std::stringstream ss;
				ss << "smallCrystal" << i << "," << j;
				G4String name = ss.str();
				G4VPhysicalVolume* smallCrystalPhys = new G4PVPlacement(0,G4ThreeVector(xPos,yPos,0),smallCrystalLog,name,crystalBlockLog,0,0);
				(void)smallCrystalPhys;
				crystalCount++;
			}
		}
	}

	// Define large crystals

	G4Box* largeCrystalSolid = new G4Box("largeCrystalSolid",2.5*cm,2.5*cm,25*cm);
	G4LogicalVolume* largeCrystalLog = new G4LogicalVolume(largeCrystalSolid,CsI,"largeCrystalLog",0,sensitiveDetector);

	int largeCrystalDimension = 36;
	for(int i = 0;i < largeCrystalDimension;i++)
	{
		for(int j = 0;j < largeCrystalDimension;j++)
		{
			double xPos = (i-largeCrystalDimension/2+.5)*5*cm;
			double yPos = (j-largeCrystalDimension/2+.5)*5*cm;
			if(fabs(xPos) > 0.6*m || fabs(yPos) > 0.6*m)
			{
				if(sqrt(pow(fabs(xPos)+2.5*cm,2)+pow(fabs(yPos)+2.5*cm,2)) < 0.95*m)
				{
					std::stringstream ss;
					ss << "largeCrystal" << i << "," << j;
					G4String name = ss.str();
					G4VPhysicalVolume* largeCrystalPhys = new G4PVPlacement(0,G4ThreeVector(xPos,yPos,0),largeCrystalLog,name,crystalBlockLog,0,0);
					(void)largeCrystalPhys;
					crystalCount++;
				}
			}
		}
	}

	G4cout << "crystal count: " << crystalCount << G4endl;

	return worldPhys;
}
